/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm.hxx"

#define UHM_ERROR_TOL 1.0e-5

int main (int argc, char **argv)
{

  // input check
  if (argc != 5) {
    printf("Try : uhmtest [n_thread][decomposition][blocksize][input_file]\n\n");
    printf(" ==== decomposition method ====\n\n");
    printf(" ** Symmetric Factorization **\n");
    printf("    UHM_CHOL      ( 1),   UHM_CHOL_PIV ( 2),   UHM_CHOL_INCPIV ( 3) \n");
    printf("    UHM_LDL_NOPIV ( 4),   UHM_LDL_PIV  ( 5),   UHM_LDL_INCPIV  ( 6) \n");
    printf("\n");
    printf(" ** General Structure Symmetric Factorization **\n");
    printf("    UHM_LU_NOPIV  ( 7),   UHM_LU_PIV   ( 8),   UHM_LU_INCPIV   ( 9) \n");
    printf("    UHM_QR        (10) \n");
    return 0;
  }

  FLA_Init();

  uhm::Mesh m;

  int n_threads, decomposition, blocksize, svd_cutoff;
  double rel_thres;
  char *filename;

  n_threads     = atoi( (argv[1]) );
  decomposition = atoi( (argv[2]) );
  blocksize     = atoi( (argv[3]) );
  filename      = argv[4];

  double t_base, t_tmp, t_build_tree, t_decompose, t_solve;
  double f_decompose, f_solve, m_estimate, m_used, m_max_used;



  printf( "BEGIN : Import mesh from file : %s\n", filename );
  m = new uhm::Mesh_;
  m->import_file( filename );
  printf( "END   : Import \n" );

  uhm::set_num_threads(n_threads);

  printf( "BEGIN : Build tree \n" );
  t_base       = uhm::timer();
  uhm::build_tree(m);
  t_build_tree = uhm::timer() - t_base;
  printf( "END   : Build tree\n" );

  printf( "New mesh info < n_nodes %d, n_elements %d >\n", 
	  m->get_n_nodes(), m->get_n_elements() );

  m->lock();

  int datatype = UHM_REAL; 
  int n_rhs    = 1;
  int is_schur = false;

  uhm::set_hier_blocksize(blocksize);

  if (decomposition < UHM_LU_NOPIV)
    m->set_symmetry(UHM_SYMMETRY);
  else
    m->set_symmetry(UHM_UNSYMMETRY);

  printf( "BEGIN : Create matrix without buffer \n" );
  m->create_matrix_without_buffer( datatype, n_rhs );
  printf( "END   : Create matrix without buffer \n" );

  printf( "BEGIN : Create buffer \n" );
  m->create_matrix_buffer(is_schur);
  printf( "END   : Create buffer \n" );

  if (decomposition < UHM_LU_NOPIV) {
    printf(" Sym Sym\n");
    m->random_spd_matrix(); 
    //m->triangularize();
  } else {
    m->random_matrix(); 
  }

  m->set_rhs();

  // ----------------------------------------------------------------
  t_decompose = 1.0e9;
  t_solve     = 1.0e9;

  printf(" ==== decomposition method ====\n\n");
  printf(" ** Symmetric Factorization **\n");
  printf("    UHM_CHOL      ( 1),   UHM_CHOL_PIV ( 2),   UHM_CHOL_INCPIV ( 3) \n");
  printf("    UHM_LDL_NOPIV ( 4),   UHM_LDL_PIV  ( 5),   UHM_LDL_INCPIV  ( 6) \n");
  printf(" ** General Structure Symmetric Factorization **\n");
  printf("    UHM_LU_NOPIV  ( 7),   UHM_LU_PIV   ( 8),   UHM_LU_INCPIV   ( 9) \n");
  printf("    UHM_QR        (10) \n");
  printf("\n\n");


  printf("BEGIN : Decomposition %d\n", decomposition);
  t_base = uhm::timer();
  uhm::decompose_with_free(m, decomposition);    
  t_tmp  = uhm::timer() - t_base;
  t_decompose = min(t_tmp, t_decompose);
  printf("END   : Decomposition\n");

  m_used      = uhm::matrix_buffer_used();
  m_max_used  = uhm::matrix_max_buffer_used();
  
  printf("BEGIN : Solve %d\n", decomposition);
  t_base = uhm::timer();
  //uhm::solve_1(m, decomposition);
  //uhm::solve_2(m, decomposition);
  uhm::solve(m, decomposition);
  t_tmp  = uhm::timer() - t_base;
  t_solve = min(t_tmp, t_solve);
  printf("END   : Solve\n");

  printf("BEGIN : Check %d\n", decomposition);
  uhm::check(m, decomposition);
  printf("END   : Check\n");

  // ----------------------------------------------------------------
  unsigned int n_dof     = m->get_n_dof();
  unsigned int n_nonzero_before = m->get_n_nonzero();
  unsigned int n_nonzero_factor = m->get_n_nonzero_factor(decomposition==UHM_CHOL);
  
  m->estimate_cost( decomposition, datatype, n_rhs,
                    f_decompose, f_solve, m_estimate );

  printf("==== Report =====\n");
  printf("Number of RHS          = %d\n", n_rhs);
  printf("Number of threads      = %d\n", uhm::get_num_threads());
  printf("Decomposition          = %d\n", decomposition);
  printf("Blocksize              = %d\n", uhm::get_hier_blocksize());
  printf("NDOF                   = %d\n", n_dof);
  printf("Non-zero entries       = input %d factored %d\n", 
         n_nonzero_before, n_nonzero_factor);
  printf("Sparsity               = %E\n", 
	 (double)n_nonzero_factor/(double)(n_dof*n_dof));
  printf("--------------------------\n");

  if (uhm::is_multithreading_enable())
    printf("Openmp is used\n");
  else
    printf("Openmp is NOT used\n");

  if (uhm::is_hier_matrix_enable()) 
    printf("Hier-Matrix is used ( blocksize = %d )\n", 
	   uhm::get_hier_blocksize());
  else
    printf("Hier-Matrix is NOT used \n");

  printf("--------------------------\n");
  printf("Time build tree (s)   = %E\n", t_build_tree);
  printf("Time tree / Time decom= %E\n", t_build_tree/t_decompose);
  printf("--------------------------\n");
  printf("Buffer estimate (MB)  = %6.0lf\n", m_estimate/1.0e6);
  printf("Buffer used     (MB)  = %6.0lf\n", m_used/1.0e6);
  printf("Max buffer used (MB)  = %6.0lf\n", m_max_used/1.0e6);
  printf("--------------------------\n");
  printf("Time decom (s)        = %E\n", t_decompose);
  printf("Time solve (s)        = %E\n", t_solve);
  printf("--------------------------\n");
  printf("FLOP decom (GFLOP)    = %6.3lf\n", f_decompose/1.0e9);
  printf("FLOP solve (GFLOP)    = %6.3lf\n", f_solve/1.0e9);
  printf("--------------------------\n");
  printf("FLOP decom (GFLOP/s)  = %6.3lf\n", f_decompose/t_decompose/1.0e9);
  printf("FLOP solve (GFLOP/s)  = %6.3lf\n", f_solve/t_solve/1.0e9);
  printf("--------------------------\n");

  if (decomposition <=10 ) {
    double residual = m->get_residual();

    printf("Residual              = %E\n", residual);
    if ( residual < UHM_ERROR_TOL ) 
      printf("TESTING UHM : **** PASS **** \n");
    else 
      printf("TESTING UHM : **** FAIL **** \n");
  }

  m->unlock();

  delete m;

  FLA_Finalize();
  return 0;
}


